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ABSTRACT 

Using images from the Hubble Space Telescope Wide-Field Camera 3, we measure the rate of diffu¬ 
sion of stars through the core of the globular cluster 47 Tucanae using a sample of young white dwarfs 
identified in these observations. This is the first direct measurement of diffusion due to gravitational 
relaxation. We find that the diffusion rate n ~ 10—13 arcsecond^ Myr“^ is consistent with theoretical 
estimates of the relaxation time in the core of 47 Tucanae of about 70 Myr. 

Subject headings: globular clusters: individual (47 Tuc) — stars: Population II, Hertzsprung-Russell 
and C-M diagrams, kinematics and dynamics 


1. INTRODUCTION 

Globular clusters have long provided an amazing lab¬ 
oratory for stellar evolution and gravitational dynamics, 
and the nearby rich cluster, 47 Tucanae, has long been 
a focus of such investigations. The key point of this in¬ 
vestigation is an interplay between these two processes. 
In particular in the core of 47 Tucanae, the timescale for 
stellar evolution and the timescale for dynamical relax¬ 
ation are similar. The relaxation time in the core of 47 


Tuc is about 70 Myr (Harris 1996). Meanwhile over a 
span of about 150 Myr the most massive stars in 47 Tu¬ 
canae evolve from a red giant star with a luminosity of 
2,000 times that of the Sun to a white dwarf with a lu¬ 
minosity less than a tenth that of the Sun. Meanwhile 
the star loses about forty percent of its mass, going from 
0.9 to 0.53 solar masses. It is these young white dwarfs 
that are the focus of this paper. 


numerous previous investigations (e.g McLaughlin et al. 

2006 

Knigge et al. 2008 

Bergbusch & Stetson|2U09 ), this 


of the Hubble Space Telescope (HST) with a mosaic that 
covers the entire core of the cluster. Probing the core of 
the cluster in the ultraviolet is advantageous in several 
ways. First the young white dwarfs are approximately as 
bright as the upper main sequence, giant and horizontal 
branch stars at 225 nm, so they are easy to find. In fact 
the brightest white dwarfs are among the brightest stars 
in the cluster and are as bright as the blue stragglers. 
Second, the point-spread function of HST is more con¬ 
centrated in the ultraviolet helping with confusion in the 
dense starfield that is the core of 47 Tuc. 

In spite of these advantages, for all but the bright¬ 
est stars, our dataset suffers from incompleteness which 
presents some unique challenges. We reliably character¬ 
ize the incompleteness as a function of position and flux 
in the two bands of interest F225W and F336W through- 
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out the colour-magnitude diagram and especially along 
the white-dwarf cooling sequence through the injection 
and recovery of about 10® artificial stars into the im¬ 
ages. Ho w we measure the completeness is described in 
detail in §3.1[ The young white dwarfs typically have a 
mass forty percent less than their progenitors, so they are 
born with less kinetic energy than their neighbors, and 
two-body interactions will typically increase the kinetic 
energy of young white dwarfs over time and change their 
spatial distribution. We introduce a simple model for 
the diffusion of the y oung white dwarfs through the core 
of the cluster ((3.2). To make the most of this unique 


dataset, we have to include stars in our sample whose 
completeness rate is well below fifty percent. We have 
developed and tested statistical techniques to character¬ 
ize the observational distribution of young white dwarfs 
in flux and space to understan d their m otion through 
the cluster and their cooling (§§3.4|3.7|) in the face of 
these potentially strong observational biases. Although 
these techniques are well known especially in gamma-ray 
astronomy, they have ne ver b een applied to stellar pop¬ 
ulations in this way, so (3.8 presents a series of Monte 
Carlo simulations to assess the potential biases of these 
techniques and verify that these techniques are indeed 
unbiased in the face of substantial incompleteness within 
the statistical uncertainties. To establish the time over 
which the white dwar fs dim we use a stellar evolution 
model outlined in §3.3[ (j^describes the best-fitting mod- 
els fo r the density and flux evolution of the white dwarfs. 
( 4.1| looks at the dynamic consequences of these results. 
( ^ outlines future directions both theoretical and obser¬ 
vational and the broader conclusions of this work. 


2. OBSERVATIONS 

A set of obs ervations with t he Advanced Camera for 
Surveys (ACS, [ Ford et al.|1998|) a nd the Wide Field Cam¬ 
era 3 (WFC3, |MacKenty| |2012p on the Hubble Space 
Telescope (HS'i'j of the core of the globular cluster 47 
Tucanae over one year provides a sensitive probe of the 
stellar populations in the core of this globular cluster 
(Cycle 12 GO-12971, PI: Richer), especially the young 
white dwarfs. Here we will focus on the observations with 
WFC3 in the UV filters, F225W and F336W. The ob¬ 
servations were performed over ten epochs from Novem- 
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ber 2012 to September 2013. Each of the exposures in 
F225W was 1080 seconds, and the exposures in F336W 
were slightly longer at 1205 seconds. Each of the over¬ 
lapping WFC3 images was registered onto the same ref¬ 
erence frame and drizzled to form a single image in 
each band from which stars were detected and character¬ 
ized, resulting in the color-magnitude diagram depicted 
in Fig. [2 

What is immediately striking in Fig. is that the dis¬ 
tribution of young white dwarfs with a median age of 
6 Myr is significantly more centrally concentrated than 
that of the older white dwarfs that have a median age 
of 127 Myr. The white-dwarf distribution appears to be¬ 
come more radially diffuse with increasing age, a signa¬ 
ture of relaxation. One concern is immediately apparent. 
The numbers of observed stars are given in legend of the 
left panel and the numbers of stars in the completeness 
corrected samples are given in the legend of the right 
panel. The sample of older white dwarfs is only about 
seventy-five percent complete on average. Furthermore, 
one would expect the completeness of these faint stars to 
be lower near the centre of the cluster, so if the complete¬ 
ness is not accounted for correctly, one could naturally 
conclude that the white dwarfs are diffusing when they 
are not in reality. In principle we would like to divide 
this sample of over 1,300 stars into subsamples some of 
which will have even smaller completeness rates. How 
can we be sure that our analysis techniques are up to the 
task of measuring this diffusion accurately in the face 
of completeness rates as low as twenty percent that vary 
dramatically with distance from the center of the cluster? 

In the following section (§ [^ we will characterize the 
completeness rate through artificial star tests, develop 
and test statistical techniques to measure the diffusion of 
white dwarfs in 47 Tucanae without binning the stars at 
all, thus preserving the maximal information content of 
these data. We will test these new algorithms on mock 
data sets that include both the completeness rate and 
flux error distribution of our sample to verify that they 
robustly determine the diffusion and flux evolution of 
the white dwarfs. The subsequent section (§13 explores 
results of these techniques on the dataset depicted in 

Fig. [3 


3. ANALYSIS 
3.1. Artificial Star Tests 

We inserted ~ 10® artificial stars into the WFC3 im¬ 
ages both in F225W and F336W over the full range 
of observed magnitudes in both bands and a range of 
distances from the center of the cluster. To determine 
the completeness rate for the white dwarfs that we have 
observed, we inserted artificial stars whose F225W and 
F336W magnitudes lie along the observed white-dwarf 
track in the CMD. The rate of recovering a star along the 
white-dwarf track of a given input magnitude in F336W 
at a given radius is the completeness rate and is depicted 
in Fig. If an artificial star along the white-dwarf track 
is detected in F336W, it is always detected in F225W as 
well. The completeness rate is both a strong function of 
radius and magnitude and is significantly different from 
unity except for the brightest stars, so accounting for 
completeness robustly is crucial in the subsequent anal¬ 
ysis. The radial bins are 100 pixels in width and the 


magnitude bins are 0.1358 wide. 

The magnitudes of the recovered stars give the error 
distribution as a function of the input magnitude and po¬ 
sition of the star in the field. Furthermore, these distri¬ 
butions are not typically normal and often asymmetric as 
well. For the analysis in §3.5| we use the cumulative distri¬ 
bution of magnitude errors as a function of position and 
input magnitude which we obtain by sorting the output 
magnitudes in a given bin and spline to obtain the cumu¬ 
lative distribution in the form of the values of the errors 
from the first to the ninety-ninth percentile. In the anal¬ 
ysis the completeness rate is interpolated over the two 
dimensions of radius and magnitude with a third-degree 
spline, and the error distributions are interpolated lin¬ 
early over the three dimensions (radius, magnitude and 
percentile). 


3.2. Diffusion and Luminosity Evolution 

Sometime during the late evolution of a turn-off star 
in 47 Tuc, the star loses about forty percent of its mass, 
going from a main-sequence star of ninety percent of a 


lar mass (Renzini & Fusi Pecci 

1988 iRenzini et al.|1996 

Moehler et al.|2UU4 Kafirai et 

af.|2UL)9). These newborn 


massive progenitors, so as they interact gravitationally 
with other stars, their velocities will increase through 
two-body relaxatio n, bringing th eir kinetic energies into 


equipartition (e.g. Spitzer 1987). Because the gravita¬ 
tional interaction is long range and the distance between 
the stars is small compared to the size of the cluster, the 
change in velocity will be dominated by distant interac¬ 
tions and small random velocity jumps, i.e. the Coulomb 
logarithm is large, ~ In where N is the number of stars 
in the core of the cluster, ~ 10®“®. These small jumps in 
velocity can be modeled as a random walk in velocity so 
the square of the velocity increases linearly in time and 

the relaxation time can be defined as t^ = [dlnu^/dt] 

In the center of the cluster, the density of stars is approx¬ 
imately constant, so the gravitational potential has the 
approximate form 




( 1 ) 


By the virial theorem the mean kinetic energy of the 
white dwarfs will equal the mean potential energy. 


\{v^) = \'^GPc{r'^) 


( 2 ) 


so the square of the distance of the white dwarfs from the 
center of the cluster will also increase linearly with time 
as a random walk; therefore, let us suppose that newly 
born white dwarfs diffuse outward through the cluster 
following the diffusion equation 


dp{r, t) 
dt 


= KV‘^p{r, t) 


( 3 ) 


where k can be related to the relaxation time as k = 
because dlnr^/dt = dhiv^/dt from Eq.This diffusion 
equation yields the Green’s function 


i(r,t) = 


8 {-KKt) 


3/2 




( 4 ) 
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Fig. 1.— Left: Color-magnitude diagram of the core of 47 Tucanae in the WFC3 filters F225W and F336W. Right: The radial distribution 
of young and old white-dwarf stars as highlighted in the color-magnitude diagram (completeness corrected numbers appear right-hand panel). 


if K is independent of time and position. This gives a 
cumulative distribution in projected radius 

C'(< i?) = (5) 

The Green’s function at t = 0 is a delta function centered 
on the center of the cluster. On the other hand, if the 
initial distribution is a Gaussian centered on the center 
of the cluster the density of white-dwarf stars near the 
center of the cluster is a function of age, t, and projected 
radius, i?, of the form 


P = p{R,t) = 


2R 


exp 




4k( t -|- to) 


4K(t + to) 

where the density distribution is normalized as 

POO 

/ p{R, t)dR — 1. 

Jo 


( 6 ) 


( 7 ) 


The dispersion of the Gaussian at t = 0 is simply given 
by 2(7^ = AkIo. Because the diffusion equation is linear, 
a sum of several Gaussians with the same value of k but 


different normalizations and values of to will also be a 
solution. 

Of course we don’t directly observe the ages of the 
white dwarfs. Rather we observe their fluxes or appar¬ 
ent magnitudes. The cooling curve of the white dwarfs 
is a relationship between time and the apparent magni¬ 
tude from the white dwarfs t(m), so the number of white 
dwarfs that we expect to observe at a given flux and 
radius is given by 

• Ot 

f{R, m) = Np{R, t{m))—C{R, m) (8) 

where N is the birthrate of the white dwarfs (assumed 
to be constant over the range of ages of the young white 
dwarfs, i.e. the past 200 Myr) and C{R,m) is the com¬ 
pleteness as a function of radius and flux. To this point 
flux errors have been neglected. 

3.3. Cooling Models 

To construct the various cooling models here, i.e. t{m) 
from Eq. we used MESA (Modules for Experiments 
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Fig. 2.— The completeness rate as function of radius from the 
center of 47 Tuc and the magnitude of the artificial star. 


in Stellar Astrophysics; Paxton et al. 2011) to perform 
simulations of stellar evolution starting with a pre-main- 
sequence model of 0.9 solar masses and a metallicity of 
Z = 3.3 X 10“^ appropriate for the cluster 47 Tucanae. 
This is sl ightly larger than the va lue for the turnoff mass 
found by Thompson et al. (2010) for the eclipsing binary 
V69 in 47 lucanae that is composed of a upper-main se¬ 
quence star of about 0.86 solar masses and a subgiant of 
0.88 solar masses. Because we are interested in the stars 
that have become young white dwarfs just recently, the 
initial masses of these stars should be slightly larger than 
the turnoff mass today. We have explicitly assumed that 
the progenitors of the white dawrfs are a uniform pop¬ 
ulation. Although there is evidence of mo dest variation 
in the chemical abundances in 47 Tuc (e.g. |Milone et al. 
2012 ), the white-dwarf c ooling sequence, at le ast at larger 


radii, appears uniform (Richer et al. 2013|. However, 
from Fig. Hit is apparerit that the core of 47 Tuc has a 
substanticd population of blue stragglers that will evolve 
to become more massive white dwarfs. Our sample has 
about 160 blue stragglers, and if we estimate the duration 
of t he main sequence for a blue straggler to be 1 — 2 Gyr 
(e.g [Sills etlTlI^OO^ , we obtain a birth-rate of blue- 
straggler white dwarfs of about 0.1 Myr“^. The number 
of giants in our field indicates a bi rth rate of about eight 
white dwarfs per million years (see Goldsbury et al.|20 12 


for further details), so the estimated contamination o: 
the white-dwarf cooling sequence is modest at about one 
percent. 

Specifically, we used SVN revision 5456 of MESA 
and started with the model lM_pre_ms_to_wd in the 
test suite. We changed the parameters initial_mass 
and initial_z of the star and adjusted the parameter 
log_L_lower_limit to —6 so the simulation would run 
well into the white dwarf cooling regime. We also re¬ 
duced the two values of the wind 77 to 0.46 (from the de¬ 
fault of 0.7) to yield a 0.53 solar mass whit e dwarf from 
the 0. 9 solar mass progenitor. Interestingly [Miglio et al.] 
(2012) argue from Kepler asteroseismic measurements of 
the stars in the metal-rich open cluster NGC 6971 that 


such values of 77 are needed to account for the mass loss 
between the red giant and red clump phases of stars in 
this metal-rich cluster. 

We defined the time of birth of the white dwarf to co¬ 
incide with the peak luminosity of the model at the tip 
of the asymptotic giant branch about 10.9 Gyr after the 
start of the simulation. This is in agreement with the 
best age of the cluster determined from main-sequence 
stars of 11.25 ± 0.21 (r andom) ±0.85 (systematic) Gyr 
Thomps on et al. ([2010 ) . Thi s age agrees with that de- 
rived by Hansen et al. (2013| from white-dwarf cooling 
(9.9±0.7 Gyr at 957o conhdence). We choose this defini¬ 
tion of the birth so that each observed white dwarf will 
have a star of similar luminosity in the cooling model. 
At this point in the evolution we have outputs from the 
MESA evolution every 100 years or so; therefore, the 
cooling curve is well sampled throughout. At each out¬ 
put time we have the value of the luminosity, radius, 
effective temperature and mass of the star. With these 
values we in terpolate the spectral models of [Trembla^^ 
et al. ( 2011 ) in surface gravity and effective tempera¬ 
ture and then scale the result to the radius of the model 


star. We use a t rue distance modulus of 13.23 (Thomp- 
son et aLl|2010|) an d a reddening of E{B — V) = 0.04 
(Safaris et al.|2007) to determine the model fluxes in the 
WFG3 band F336W. We used the standard extinction 


curve of Fitzpatrick (1999) with Av/E{B — V) = 3.1. 
We have purposefully used a distance and reddening de¬ 
termined from main sequence stars to avoid a potential 
circularity in usi ng the white dwarf mo dels themselves to 
fix the distance. Woodley et al. ([2012 ) inferred a slightly 
larger true distance modulus of 13.36 ± 0.02 ± 0.06 from 
the white-dwarf spectral energy distributions. 

The brightest white dwarf in our sample has F336W = 
14.92. According to the models this corresponds to 
an age of 110,000 years, an effective temperature of 
100,000 K, a luminosity of 1,600 Lq, and a radius of 
0.13 Rq. Its mass is 0.53 solar masses. The faintest 
white dwarf in our sample has F336W = 25.4, yielding 
an age 1.2 Gyr, an effective temperature 8,700 K, a lu¬ 
minosity of 10“^ Lq and a radius of 0.013 R©, one tenth 
of the radius of the brightest white dwarf. Clearly the 
brightest white dwarf in our sample is not a white dwarf 
in the usual sense because thermal energy plays an im¬ 
portant role in the pressure balance of the star. For this 
brightest star log <7 = 5.93 which is less than the min¬ 
imum of the atmosphere model grid (log g = 6) so we 
have to extrapolate slightly off of the gri d, bu t only for 
this brightest star. For the simulations in §3.8| we did not 
use this particular model, but similar ones of the same 
white dwarf mass with different neutrino cooling rates or 
initial metallicities also generated with MESA. 

3.4. Likelihood Function 

The model outlined in §3.2| predicts the number of 
white dwarfs as a function of magnitude and position. 
Let us divide the space of position and magnitude into 
bins of width AR and Am and where the bins are num¬ 
bered with indices j and k respectively. The probability 
of finding n stars in a particular bin is given by 


P{n; f{Rj,mk)) = 


[/(Rj, TOfc)ARATO]" 


n\ 


(9) 
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Now we imagine dividing the sample into so many bins 
that there is either a single star in a bin or no stars at 
all, we have 

J/(Pj,mfe)APAm one star 
I 1 no star ^ ' 


We can define the likelihood as the logarithm of the prod¬ 
uct of the probabilities of observing the number of stars 
in each bin. Since the bins are so small we can replace Rj 
and rrik for the bins with stars in them with the measured 
values for that particular star Ri and rrii. This gives the 
so-called “ unbinned likelihood” of observing the samp le 
as follows ( |Cas h| |1979[ [Mattox et al.||199^ |Davis|[2M4| ) 


logL = ^log/(i?i,mi) -'^f{Rj,mk)ARAm. (11) 


j,k 


We have dropped the constant widths of the bins from 
the first term which is a sum over the observed stars; 
consequently, the absolute value of the likelihood is not 
important, just differences matter. The second term is a 
sum over the really narrow (and arbitrary) bins that we 
have defined, so we have 


'^f{Rj,mk)ARAm = j f{Rj,mk)dRdm = Np^^d 

j,k 

( 12 ) 

where Np-^^^ is the number of stars that the model pre¬ 
dicts that we will observe, so finally we have 

logL = ^log/(Pi,mj) - Apred (13) 

i 

where the summation is over the observed stars. The 
integral for Np^ed when combined with Eq. yields 

pti ^Tmax 

Npred = N / p{R,t)C{R,m)dRdt (14) 
Jo Jo 


3.5. Magnitude Errors 

An important complication to the analysis is that the 
measured magnitudes are not the same as the actual 
magnitudes of the stars; in particular the error distribu¬ 
tion is not normal or even symmetric. This transforms 
the model distribution function via a convolution. 


/ OO 

f{R,m')g{R,m',m — m')dm' (18) 

-OO 

/ OO 

f{R, m — Ain)g{R, m — Am, Am)d (Am) 

-OO 

= f f{R,m — Am)dG (19) 

Jo 


where 


/ Am 

g{R,m—Am, Am')d{Am') (20) 

-OO 


is the cumulative distribution of magnitude errors with 
the observed radius and magnitude fixed. If we calcu¬ 
late the percentiles of the magnitude error distribution 
as Amj we can approximate the integral as the sum 

/'(P,to) = ^^/(P,TO-Atoj), (21) 

j 


SO for a given star i we have 


f'{Ri,mi)=N'^p{R,t{mi - Amj)) x 
j 

t' {mi — Amj) C{Ri,mi — Amj) (22) 

where t' = dt/dm. Thi s new function f'{Ri,mi) can be 
substituted into Eq.f^to yield a likelihood including the 
magnitude errors. We will assume that the magnitude 
errors do not affect our estimate of Ap^ed; this simplifies 
the analysis. We will verify our technique with Monte 
Carlo simulations in §3.8| 


or 


Apred = A 


' mo 'J 0 


p{R,t)C{R,m)-^dRdm. (15) 
dm 


If we take the luminosity function as fixed and try to 
maximize the likelihood with respect to the diffusion 
model 


logT = ^log 


Np{R,t{mi)) 


= ^log Np{R,t{mi)) 

i 

dt 


dm 




G {Ri , rjii) 


-A 


^log 


dm 


C{Ri,mi) 


-A, 


pred 


pred 


(16) 


(17) 


where the second summation does not depend on the dif¬ 
fusion model so it is constant with respect to changes in 
the diffusion model and can be dropped from the loga¬ 
rithm of the likelihood. However, it must be included if 
one wants to compare different cooling curves, t{m). 


3.6. Constraining the luminosity function 

We can construct a maximum likelihood estimator of 
the luminosity function of the white dwarfs or alterna¬ 
tively the cooling curve as follows 

= Np{R,t)^^ AiS{m — mi) (24) 


where i is an index that runs over the observed stars. 
With this model we can define a likelihood function for 
the stars that we observe 


log L = ^ log NAip{Ri, ti) 


-A 


pred 


(25) 


where a multiplicative constant (infinite in this case) and 
the completeness for each star have been dropped from 
the logarithm. 

Substituting the trial luminosity function Eq. |23| yields 

r^max 

Apred = VA,A/ p{R,ti)C{R,mi)dR. (26) 
, Jo 



















6 


If we maximize the likelihood with respect to the values 
of Ak we obtain 


d\og[p{Ri,U)C{Ri,mi)\ d log [p{Ri,ti)] dU 


dAk 


dU dAk 


where 


so 


dtj 

dAk 


1 if k <i 
0 if k > i 


(27) 


(28) 


d\og[Aip{Ri,ti)C{R^,mi)] _ 1 91og [p(i?i, t^)] 

^ dAk ~ Ak ^ ' 


and 


9 log lp(R,t)] 
dt 


i>k 


dti 


R^ 


Ault + tg)^ ^ "b ^0 


(29) 

(30) 


Taking the derivative of Np^ed yields the second part of 
the variance in log L, 


f^JSJ p'^max 

= N / p{R, tk)C(R, mk)dR + 
Jo 

“ dp{R,ti) dU 

dti dAk 


Y^a^n 


(31) 

C{R, mi)dR 


= N p{R,tk)C{R,mk)dR + 


(32) 


Ya.n 


i>k 


dp{R, tj) 

dti 


C{R, mi)dR. 


Combining these results with dlogL/dAk = 0 yields an 
equation of the form 


1 ^ P'^max 

— =N p{R,tk)C{R,mk)dR + 

Ak Jo 


(33) 


E 

i>k 


a,n 


dp{R,U) 


dti 


C{R, mi)dR — 


a log [p{R^Ai)\ 
dU 


or a matrix equation of the form 


d^ikAi — bk + 

J±h 


where 


and 


10 11 i < k 


(34) 


(35) 


M,=N 


dp{R,ti) 


dti 


C{R,m^)dR. (36) 


The vector bk is given by 


^>k -^0 

(37) 


Although this matrix equation has as many rows as there 
are stars in the sample, it is straightforward to solve at 
least formally in two ways. The first is 


Ak — 


1 


AIikAi — bk 


(38) 


and the second is 


A,; = 


bi + {Ai) ^ bi+i + (Ai-|_i) 


M, 


M, 


, An — 


(A.)- 


i+l 


Mn 


(39) 

The values of bi and Mi, of course depend on the values 
of Ai through the parameter tk , so the solution must pro¬ 
ceed iteratively perhaps while minimizing with respect to 
the other parameters of the model k and tg- 
For each value of the diffusion parameters, we chose to 
iterate Eq. three times to determine the values of Ak 
with in a loop of two iterations where Mi (Eq. 361 and bk 


(Eq. 371 vary. Given this new trial luminosity tunction 
the emulsion parameters are varied to find the maximum 
likelihood, and the iterative solution of the luminosity 
function is repeated. These two steps are repeated until 
the values of the diffusion parameters from one iteration 
to the next have changed by less than one part per hun¬ 
dred. 

An interesting limit is when the density distribution 
is independent of time. This understandably yields a 
simpler solution for A^. In particular, Mi = 0 so 


A^ = -^ = 

bi 


N / p{R,mi)C{R,mi)dR 


(40) 


where the underlying density distribution is normalized. 
The weight is not the reciprocal of the completeness for 
star i but rather the reciprocal of the mean of the com¬ 
pleteness of a star with the flux of star i over the density 
distribution. The latter could be evaluated by taking 
the mean of the completeness measured for all the stars 
in the sample in a magnitude range about star i suf¬ 
ficiently wide to sample the density distribution. It is 
important to note that the weight is the reciprocal of 
the mean of the completeness not the mean of the recip¬ 
rocal. If the completeness does not depend strongly on 
radius, these two will approximately coincide. Finally 
if the density distribution is not known a priori and is 
not modeled, the weight for a particular star is simply 
given by Ai = \C{Ri,mi)Y. We call this “Inv Comp” 
in Figs. and [l^ 

The likelihood is invariant under changes in the birth 
rate of the white dwarfs {N) if one also changes the values 
of K, Ai and to as follows: 


N —>■ aN, Ai 


, -1 


Ai,to —t a ^tg and k — an. (41) 

That is the time scale cannot be fixed without some ad¬ 
ditional input such as a theoretical cooling curve or an 
independent estimate of the white dwarf birthrate. The 
quantities AiN, k/N and Nto are invariant with respect 
to this transformation. In our dataset when we use this 
modeling technique, we fix the value of N to the value 
inferred by the nu mber of giants in our field as in [Golds-] 
bury et al. (2012). 


3.7. Constraining the luminosity function with errors 
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We s tart the analysis including magnitude errors with 
Eq. [island which when combined yield, 

/ OO _ 

Np{R, t)C{R, m') ^ Ai5{m' — rm) x 

i 

g{R,m',m — m')dm' (42) 

= N'^Aip{R,ti)g{R,mt,m- nil) (43) 

i 

With this model we can define a likelihood function for 
the stars that we observe 


logE = ^log 


N'^Aip{RjAi)giRj,nii,mj - mi)C{Rj, 


-iV, 


pred • 



Note how the magnitude error essentially translates into 
a spread in the age of the observed stars. 


9log7^ _ ■sp I p{Rj,ti)9{Rj,nii,mj - mi)C{Rj,mi) 

dAk ^:^f^\j2i^iPiRjAi)9{Rj,nii,mj - mi)C {Rj, nii) 


X 


5ik + Ai 


9log [p{Ri,U)] 
dti 


^A^pred 

dAk 


(45) 


Although we have included this additional complication 
in the derivations for completeness, we have found that 
the inclusion of error convolution in modeling simulated 
data does not affect the fitting results, so we did not 
include this in the modeling of the dynamics while si¬ 
multaneously determining the luminosity functions. 


3.8. Monte Carlo Simulations 



To test these techniques in the face of the challenges of 
incompleteness and magnitude errors present in our data, 
we simulated typically on the order of 10,000 catalogs of 
the same size as our dataset with a known luminosity 
function and a known diffusion model and attempted to 
recover the input parameters. In both cases, the age of 
the star is selected first to be between zero and 1.5 Gyr. 
Given this age the model cooling curve determines the 
F336W magnitude. Second, a radius is selected from 
the cumulative distribution in projected radius (Eq®. 
Given the radius and magnitude of the candidate for the 
catalog, the completeness for this star is calculated and 
the star is included in the sample with this probabil¬ 
ity. Finally, the magnitude errors are applied by draw¬ 
ing from the magnitude error distribution. We created 
a sample of 3,167 stars — the same as in the WFC3 
white-dwarf sample. The fitting procedure followed two 
different strategies. 

The first was to assume a fixed cooling curve and try 
to find the density evolution to determine whether the 
process is biased in determining the diffusion parameters 
and the typical errors. Finally, we performed simulations 
where we did not convolve the models with the error dis¬ 
tribution to calculate the likelihood (in all cases errors 
were applied to the simulated data) to see whether the 
omission of this step introduced biases. The second strat¬ 
egy did not assume a cooling curve and determined the 
cooling curve as a part of the process of determining the 
diffusion. We did not convolve the cooling curve with the 


Fig. 3.— Upper panel: The values of the fitted parameters are 
typically unbiased with respect to the input values in the simu¬ 
lations, here N is depicted. The input values for the two types 
of simulations are given by the vertical lines. Lower panel: The 
values of tg. 


TABLE 1 

Parameters from the Monte Carlo Simulations: Means, 
Standard Derivations, Input Values 


Technique 

/€ 

Input 

to 

Input 

Full Modeling 

3.4 ±0.3 

3.58 

560 ± 80 

531 

No-Error Convolution 

3.7 ± 0.4 

3.71 

530 ± 80 

515 

Unfixed LF 

7.5 ± 0.6 

7.26 

220 ± 30 

231 



N 

Input 

logL 

Full Modeling 

5.42 ± 0.03 

5.44 

±50 

No-Error Convolution 

5.45 ± 0.03 

5.45 

±22 

Unfixed LF 

- 

±20 



error distribution while fitting the model; however, the 
fake catalogs were created in the same way as in the first 
strategy. In this technique the resulting co olin g curve 
can be multiplied by a constant factor (Eq. 41), so we 
determine the values of k and tg by fixing the value of N 
to the one used to build the catalog. This also fixes the 
age estimates of all of the white dwarfs in the sample. 

The results of these simulations are depicted in 
Fig. [3] and 1^ and in Tab. The key results of the sim- 










































Fig. 4.— Upper panel: The values of the fitted parameters are 
typically unbiased with respect to the input values in the simu¬ 
lations. Here k, is depicted. The input values for the two types 
of simulations are given by the vertical lines. Lower panel: The 
distribution of log L is significantly wider when error convolution 
is included in the fitting process. 


Fig. 5.— Upper panel: Monte Carlo simulations of an evolving 
density distribution that are fit by an evolving distribution yield 
an unbiased estimate of the cooling curve. Lower panel: If the 
evolving density distribution is not included in the fitting the re¬ 
sulting cooling curve is typically steeper than the input cooling 
curve; however, this difference is subtle. The initial guess is the 
starting point for the iterative solution of the cooling curve. 


ulations are that the likelihood fitting of the diffusion 
model results in an unbiased estimate of the diffusion 
parameters regardless of whet her t he fitting technique 
includes the magnitude errors (j |3.5[ ). Furthermore, even 
when one fits for the luminosity function as well one can 
obtain reliable estimates of the diffusion model without 
prior knowledge of the cooling curve; of course, in this 
latter case the timescales of the diffusion rely on an in¬ 
dependent estimate of the birth rate of the white dwarfs 
N. Observationally, this is determined from a sample of 
giant stars numbering in the thousands (see Fig. y so 
the statistical error in this determination is small. Uyp- 
ically the birth rate is recovered with an uncertainty of 
less than one percent and the diffusion rate with an un¬ 
certainty of ten percent and to with an uncertainty of 
fifteen percent. The errors in to and k are correlated so 
the error in toK is typically less than ten percent. 

In the second type of simulation, we found the density 
evolution along with an estimate of the cooling curve, 
so this cooling curve can be compared with the input 
cooling curve for the simulations. Furthermore, the de¬ 
termination of the cooling curve is iterative, so we have 


to give an initial guess of the curve. The input, the ini¬ 
tial guess and the results are given in Fig.[^ We can also 
fit for just the cooling curve and assume that the den¬ 
sity distribution does not evolve or not assume a density 
model a t all and use the per star completeness as out¬ 
lined in §3.6| Fig. [^highlights the difference between the 
model age and the inferred age with the various likeli¬ 
hood techniques. For young white dwarfs the uncertain¬ 
ties are large (because there are few young white dwarfs 
in the sample), but for old white dwarfs there is a small 
bias of order of ten percent in the inferred age, the sign 
of which depends on the technique. Again this is on the 
order of the relative errors in the diffusion parameters. 


4. RESULTS 

The results of the diffusion model fitting are given in 
Tab. The results do not depend strongly on the mod¬ 
eling technique, especially the assumed cooling curve for 
the white dwarfs. The inferred relaxation times are also 
in goo d agreement with the value tabulated by |Harris| 
(|1996|). Fig. [t] shows the posterior probability distnEm 
tion tor the various parameters and how the uncertainties 
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Fig. 6.— The relative differences between the first, second and 
third quartiles and the input model for the different luminosity 
function fitting techniques. 


TABLE 2 

Diffusion parameters from the likelihood fitting. The 

VALUES OF (TO AND tr ARE GIVEN BY ^2K,tQ AND r^/K WHERE 

rc = 22". The posterior probabilities for the full model 

ARE GIVEN IN FiG. [3 


Model 

[(//)2 Myr-1] 

*0 

[Myr] 

^^0 

n 

N 

[Myr-l] 

tri'f'c) 

[Myr] 

Full 

13.1 

166 

66 

7.07 

37 

No Errors 

13.1 

166 

66 

7.07 

37 

2 Gaussians 

12.8 

14.9 

19.5 

1.90 

38 

(No Errors) 


260 

82 

5.32 


Free LF 

9.80 

241 

69 

8.10 

49 


are correlated with each other. An important conclusion 
is that the no-diffusion model [i.e. tt = 0 ) is excluded at 
high confidence. 

To find whether we could better fit the radial distribu¬ 
tion with a sum of Gaussians, we performed the fitting 
with two and three Gaussians. We did not include the 
error convolution in the fitting models. The fit with two 
Gaussians has a value of logL that is lower by 75 from 
a fit with a single Gaussian. From Fig we can see 
that this is a significantly better fit. However, the de¬ 
crease in logL by adding a third Gaussian is only 2; fur¬ 
thermore, the third Gaussian has a very low value of N 
so it does not affect the resulting distributions strongly. 
Tab.j^shows that the diffusion parameters from the two- 
Gaussian fit only differ slightly from the one-Gaussian 
fits. In any case these differences lie within the sta¬ 
tistical errors. We can compare the best-fitting model 
density distributions as a function of time with the ob¬ 
served (completeness corrected) density distributions for 
several age ranges of white dwarf. The diffusion model 
for the median age of the white dwarfs in each bin is de¬ 
picted with a solid line for the one-Gaussian model and 
a dot-dashed line for the two-Gaussian model. The two- 
Gaussian model does a better job at following the dis¬ 
tribution of the white dwarfs especially at smaller radii. 


4.1. Two-Body Relaxation 


Fig.H] depicts the radial distribution of white dwarfs 
of various ages. Each bin is 50 Myr wide, and the bins 
are centered on 25 Myr, 125 Myr and 225 Myr. The 
evolution at up to a few core radii (about 60 arcseconds) 
is dramatic from 25 to 125 Myr and modest thereafter. 
Outside 60 arcseconds the cumulative distributions are 
nearly parallel indicating little evolution in this region 
at early times. The simple diffusion models used here 
assume that the diffusion coefficient is constant in space 
and in time, so the models continue to evolve at late time 
and for all radii. At the smaller radii the white dwarfs 
reach the distribution corresponding to their masses after 
about 100 Myr and stop diffusing. 

Fig. 1^ focuses on the outer half of the WFG3 field. 
Here we see more evolution between the second and 
third epochs with little early evolution. This indicates 
the increase in the relaxation time as the stellar density 
decreases. The white dwarfs diffuse modestly over the 
first 100 Myr and more dramatically during the second 
200 Myr. The white dwarfs as expected from theoretical 
considerations suffer diffusion that is a function of radius 
and time and beyond the scope of the simple model used 
to quantify the diffusion in this paper. However, this 
model does capture the diffusion within a few core radii 
for a few core relaxation times. 


5. CONCLUSIONS 
5.1. Further analysis 

In this paper we used the Green’s function (Eq. to 
model the diffusion of the stars through the cluster. We 
simply took the initial conditions to be a Gaussian or a 
sum of Gaussians centered on the center of the cluster. 
This allowed for a simple closed-form expression for the 
density function in spherical coordinates and in projec¬ 
tion as well. Without relaxing the spherical symmetry 
one could imagine much more general initial conditions. 
In fact we have an estimate of the initial conditions in 
the form of the projected radial distribution of the stars 
on the upper main sequence. This distribution could be 
possibly deproje cted as a lowered-isotherm al distribution 
in phase space (|Michie||1963| |King||1966|) and convolved 
with the Gaussian Green’s function, Fq. to give the 
expected density distribution as a function of time. This 
technique shares the advantage of the technique used in 
this paper that the density distribution can be guaran¬ 
teed to be positive because the convolution of the positive 
kernel with a positive distribution is necessarily positive; 
however, the density distribution even in spherical coor¬ 
dinates is not available in closed form. 

A second strategy would be to expand the initial den¬ 
sity distribution in terms of spherical Bessel functions 
and spherical harmonics. If we restrict ourselves to an 
initially spherical distribution we have 


p{r,t)= f dka{k)e-'^^^^^-^^^ (46) 

Jo 

where the coefficients a(fc) are determined from the initial 
density distribution 

a{k) = - r dr {krf p{r, 0 ) . ( 47 ) 

TT Jq kr 

If the initial density distribution can be well represented 
with a few values of k, then the density evolution is 
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Fig. 7.— The posterior probabilities of the best-fitting parameters. The upper panels depict the covariance among the various parameters. 
The contours trace probabilities a factor of e, and smaller than the maximum likelihood. The color scale also gives the natural logarithm 
of the relative likelihood. The lower panels give the likelihood as a function of a single parameter with the other parameters integrated out. 



Projected Radius [arcseconds] 


Fig. 8 .— Radial distribution of white-dwarf samples of various 
age ranges and median ages (given in parenthesis) and the best¬ 
fitting tw o-Gaussian di ffusion models superimposed. The core ra¬ 
dius from [Harris] | [T99^ is depicted by the vertical line. 

straightforward to evolve forward and backward in time; 
however, it is no longer guaranteed to be positive even 
at the initial time if only a range of values of k are con¬ 
sidered in a{k). 

From the point of view of the likelihood analysis, a 
natural next step would be to use the additional infor¬ 
mation available with the current observations, i.e. the 
flux in the F225W band. This would provide an addi¬ 
tional constraint on the ages of the white dwarf stars or 
alternatively constrain the cooling curve in both bands. 
In the first case one would perhaps get better constraints 
on the dynamical evolution and could also fit for the dis¬ 
tance and reddening to the cluster and possibly the mass 
of the white dwarfs or specifics of the cooling mechanism. 



Fig. 9.— Radial distribution of white-dwarf samples of various 
age ranges and median ages (given in parenthesis) and the best¬ 
fitting two-Gaussian diffusion models superimposed for the outer 
half of the region. 

In the second case one would get a cooling curve in a sec¬ 
ond band. It is straightforward to see that the weights 
for the cooling curve in F225W would be the same as 
in F336W, so simply plotting the inferred ages of the 
white dwarfs from Fig. [T0| against the F225W magnitude 
would yield the cooling curve in F225W. The agreement 
with the F336W model is poorer at early times but im¬ 
proves with age and lasts until nearly I Gyr. In the 
context of this paper, we obtain similar diffusion param¬ 
eters whether we fit a luminosity function or assume a 
theoretical model. 

5.2. Theoretical directions 
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Fig. 10.— Flux in the WFC3 band F336W as a function of 
time since the peak luminosity of the star that we define to be 
the birth of the ‘white d warf. The model curve assumes a true 
distance modulus of 13.23 jThompson et al.|2010b and a reddening 
oi E{B—V) = 0.04 (|Salaris et al.|2Uu71l. 4 he *'lnv Comp” technique 
ignores the effects ot dittusion in modeling the stars and uses the 
completeness rate corresponding to the observed magnitude and 
radius of each star. 


As argued in §3.2| the interactions with other stars 
cause the white dwarfs to diffuse in velocity not in ra¬ 
dius. However, we argued using the virial theorem that 
this diffusion in velocity would be manifest as a diffusion 
in radius as well. Furthermore, our simple model as¬ 
sumes that the diffusion coefficient is constant in space 
and time when in fact with time the white dwarf distri¬ 
bution approaches that of stars of similar mass so the 
diffusion must cease and also at larger radii the diffusion 
must happen more slowly. We see both of these effects in 
Fig. 1^ and 1^ How this diffusion actually manifests itself 
coula be simulated in two possible ways. 

The first is direct numerical simulation of on order of 
one million stars that form the central regions of the glob¬ 
ular cluster 47 Tucanae. Although on the face of it, this 
appears to be a Herculean labor when the state of the art 
direct calculation of the two body interactions in a globu¬ 
lar cluster involve merely ^10® stars and the simulation 
in question would normally take 100 times longer. How¬ 
ever, we are only interested in the dynamical evolution of 
the young white dwarfs over about one hundredth of the 
age of the cluster (100-150 Myr out of 10 Gyr). Secondly, 
because we are not interested in the long term evolution 
of the cluster, neither stellar evolution nor the dynamics 
of binaries should play an important role in this process. 


Bergbusch, P. A. & Stetson, P. B. 2009, AJ, 138, 1455 
Cash, W. 1979, ApJ, 228, 939 

Davis, D. 2014, Likelihood Tutorial, http://fermi.gsfc .nasa. 
gov/ ssc/data/6uialysis/scitools/likeiihoocl_tutorial .html 
accessed: 2U14-U8-11 
Fitzpatrick, E. L. 1999, PASP, 111, 63 


These two simplifications result in a factor of a thou- 
sand speed up to obtain results and these calculations 
are already underway. 


The second direction would be to model the diffusion 
in phase space using a Fokker-Planck or Monte Carlo 


scheme (IGiersz h Heggie|2011 

Pattabiraman et al. 12013 

Hong et ai.|2U13 

). yuch simulations would be more rapid 


than a direct n-body simulations and possibly yield more 
physical insights. 


5.3. Further observations 

Fo llowing the arguments of the preceding subsection 
§5.2| a natural direction would be to measure the proper 
motions of the white dwarfs in the core of 47 Tuc with a 
second epoch of observations. Because we already have 
the colors of the white dwarfs, only observations in a 
single band would be required and possibly not as deep 
as the present set of observations because the stars have 
already been detected. To obtain the most precise po¬ 
sitions and to minimize the crowding, the bluest band 
would be best, i.e. F225W, and possibly over only a 
portion of the field of the current data, because here the 
goal would be to verify the current result by finding the 
corresponding signal in velocity space, so a full sample 
of 3,000 plus white dwarfs may not be required. 

5.4. Final remarks 

We have measured directly for the first time the dy¬ 
namical relaxation of stars in a globular cluster. To do 
this we have introduced new statistical techniques for 
the characterization of stellar populations. These tech¬ 
niques can robustly and straightforwardly account for 
high incompleteness and non-Gaussian magnitude errors. 
They can be applied to a wide variety of questions from 
globular cluster dynamics to galaxy luminosity functions. 
There are many avenues for further investigation such as 
a more thorough analysis of the existing data using the 
information from the second band, the simulation of the 
relaxation of young white dwarfs in numerical models and 
measuring the proper motions of the young white dwarfs 
to search for signatures of relaxation in their velocities 
as well. 

This research is based on NASA/ESA Hubble Space 
Telescope observations obtained at the Space Telescope 
Science Institute, which is operated by the Association 
of Universities for Research in Astronomy Inc. under 
NASA contract NAS5-26555. These observations are 
associated with proposal GO-12971 (PI: Richer). This 
work was supported by NASA/HST grant GO-12971, 
the Natural Sciences and Engineering Research Gouncil 
of Ganada, the Canadian Foundation for Innovation, the 
British Columbia Knowledge Development Fund. This 
project was supported by the National Science Founda¬ 
tion (NSF) through grant AST-1211719. It has made 
used of the NASA ADS and arXiv.org. 
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